The Shape and Stabihty of a Viscous Thread. 
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When a viscous fluid, like oil or syrup, streams from a small orifice and falls freely under gravity, 
it forms a long slender thread, which can be maintained in a stable, stationary state with lengths up 
to several meters. We shall discuss the shape of such liquid threads and their surprising stability. 
It turns out that the strong advection of the falling fluid can almost outrun the Rayleigh-Plateau 
instability. Even for a very viscous fluid like sirup or silicone oil, the asymptotic shape and stability 
is independent of viscosity and small perturbations grow with time as exp(C t^), where the constant 
is independent of viscosity. The corresponding spatial growth has the form exp((z/L)8), where z is 
the down stream distance and L ~ Q^a~^g and where a is the surface tension, g is the gravity and 
Q is the flux. However, the value of viscosity determines the break-up length of a thread Li, ~ u^'''^ 
and thus the possibility of observing the exp{Ct^) type asymptotics. 

PACS numbers: 47.20.-k,47.20.Cq,47.20.Dr,47.20.Gv,47.15.-x,47.20.Ky,47.54.+r 



I. INTRODUCTION 



When honey or sirup is poured from an outlet, one 
easily generates very long threads of flowing fluid of sur- 
prising beauty and stabihty. A uniform column of fluid 
is unstable due to surface tension effects - the famous 
Rayleigh-Plateau instability . Viscosity diminishes the 
strength, but does not remove the instability, and thus 
the observation of stable falling viscous threads of, say, 
two meters is surprising. In the present paper we shall 
discuss the shape and stability of such falling viscous jets 
or threads. We should note from the outset that we are 
confining our attention to Newtonian fluids (e.g. sirup or 
Silicone oil). 

Our starting point is the lubrication approximation 
(see. e.g. j^), which only takes into account the leading 
order dependence of the velocity field on the radial vari- 
able and of which we give a short derivation in section 
II. In section III we study the stationary solutions, and 
in particular their asymptotic forms. The final asymp- 
totics (for large downstream distance z) is always gov- 
erned solely by gravity, as in a free fall. Then we pro- 
ceed with linear stability analysis (section IV- IV). After 
a recapitulation of the classical Rayleigh-Plateau insta- 
bility in the lubrication approximation in the absence of 
gravity, we then study the full linear stability problem 
of a falling thread using an Eulerian description in the 
comoving frame. We consider first the inviscid regime, 
which determines the asymptotics for large times if it ex- 
ist at all. We then asses the importance of viscosity and 
find that it determines the opposite asymptotics of spa- 
tially growing modes at short times. Small viscosity thus 
means that the flow breaks up without ever reaching the 
inviscid asymptotic state. 
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II. DERIVATION OF THE MODEL. 

To derive our model we use the lubrication approxi- 
mation, and our derivation is thus very close to the one 
given e.g by Eggers Since we wish to include gravity 
our equations are scaled differently, however. The veloc- 
ity field is assumed to be axisymmetrical, and it is con- 
venient to use cylindrical coordinates. We thus assume 
that the velocity field has the form: v = ucr + wez, where 
r is the radial coordinate and z is the vertical coordinate, 
measured positive downwards. We assume that the veloc- 
ity field has no azimuthal component. The Navier-Stokes 
equation and continuity equation 

Ut + UUr + WUz = —Pr/p+t^{{rUr)r/r + Uzz—u/r'^) 
Wt+UWr+WWz = —pr/p + g + l'{{rWr)r/r + Wzz){l) 
(ru)r + Wz — 

and we assume that the fluid is conflned to a thin axi- 
ally symmetric thread with a free surface at r = h{z,t). 
We expand pressure and velocity fields in power series 
in r, and assume that the expansion parameter r is 
(asymptotically) small with respect to vertical coordi- 
nate z at a given cross-section of the fluid thread, i.e. 
r/z — > 0, z ^ oo: 

w = wo{z,t) + W2{z,t)r'^ + . . . 

u = -~woz{z,t)r/2-W2z{z,ty + ... (2) 

p = pQ{z,t)+p2{z,t)r'^ + ... 

Here the expression for u guarantees that the velocity 
field is divergence-free. 

Inserting this expansion into the Navier-Stokes equation 
gives to leading order 

Wot + wqWqz = -Poz/p + 9 + v{wozz + 4^2) (3) 

To close the equation, we need to express W2 and po in 
terms of wo by using the dynamic boundary condition. 

an — —ann (4) 



2 



where ct is a stress tensor, k is a mean curvature of the 
surface, a is the coefficient of surface tension and n is a 
unit normal vector, pointing into the fluid. In terms of 
the fluid surface r — h{z,t), the normal vector is n = 

{nr,n,) = {-l,h^)/y/T+M- 

The only nonzero components of stress tensor are (see 
e.g. i): 

cTzz = -P + 2iypwz = -pa - vpiVQz + ••• (5) 
arz = vp{u^+Wr)^{2w2~WQ^j2)r+ ... 

and inserting into Q 



arrflr + (Jrznz — —aKTlr 
azrlT-r + CTzznz = —aKUz 



(6) 



gives (canceling the common multiplier 1/^1 + /i^) 



Po + vpwoz + vpi'2w2 - ■WQzz/'2)hh^ = aK (7) 
vp{2w2 - WQzz/2)h ~ {-po + 2i'pwaz)hz = anh^ (8) 

Neglecting, again for a thin thread, the nonlinear term 
hhz we obtain 

Po = - vpwoz (9) 

3 1 

W2 = -jWozrz/r + ^wq^^ (10) 

Using these expressions for po a-nd ?i'2 in (|2Jl we get: 

Wot +«^oW'02 = ~~'^0z +g + v ^Swozz + 6woz-^^ (11) 

The kinematic boundary condition leads to the conser- 
vation law for the cross-section of the thread: 



ih'^ + {woh^)z = 0, 



The curvature is: 



(12) 



(13) 



but since we are here interested in asymptotic properties 
of thin threads, we neglect the curvature in the (r, z) 
plane compared to the one around the axis of the thread, 
and assume that /i^ <C 1. Thus we shall, throughout the 
paper use the approximation: 



(14) 



Thus the our model has the following form (using a 
a/p): 



{h'')t + {woh^)z = 0, 



(15) 



wot+WQWoz = ^^1^) +5 + 3i^ (16) 



We now introduce dimensionless variables through 

z az, t ^ (3t /i — !■ ah, wq — > —v 

P 



(17) 



where a and (3 are dimensional coefficients. Thus (|16|) 
acquires the following form: 



/32 /32 fi 
vt + vvz = — g ^cr - 



3t^-T7 , o (18) 



whereas pSf) preserves its form since it is homogeneous 
in space and time variables. 

We choose a and fi such that the two first coefficients 
on the RHS of is unity, i.e. : 



1 _i 
a = a'^ g ^ 



j3 — a^g 4 



(19) 



This allows us to consider both viscid and inviscid 
cases by means of the last coefficient 7 = "ivfijo? = 
'iva~^^^g~^/^ . Note that this choice of rescaling means 
that lengths are measured in units of the capillary length 
If, — a2g^2 — a. Thus the non-dimensionlized model 
has the following form: 



{h% + ivh^)z = 



Vt + VVz 



(20) 



Typical values for 7 are jsirup ~ 100 (with simi- 
lar values for heavy silicone oils), ^glycerol ~ 0.4 and 
« 0.004. 
Solution 



III. STATIONARY SOLUTIONS. 

The shape of a stationary thread has been studied 
by several authors (see |^]-^]), but since the results are 
somewhat scattered and incomplete, we have found it im- 
portant describe the stationary states in some detail. For 
stationary solutions the non-dimensional flux: q = h^v is 
constant and we end up with the following equation for 
velocity field only: 



1 



qv 



1- 



(22) 



It is possible to remove dependence on q by an appro- 
priate rescaling, but we prefer to keep it since the total 
flux is the only parameter of the model that can easily 
be changed in a typical experiment. The flux q is by the 
scaling ()19(l related to the physical flux Q as 



(23) 



made 

When H22(l is solved forward in z, i.e. as an "initial 
value problem" , the typical solutions will diverge for large 
z. This can be circumvented by integrating backwards 
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noting that the fixed point {v,Vz) = (0,0) has a well- 
defined unstable manifold (separating solutions that di- 
verge to plus or minus infinity), which upon backward 
integration becomes a stable manifold. In Fig. ^ we 
show typical phase space trajectories found by solving 
()22|l numerically by means of a fourth-order Runge-Kutta 
method, starting from "initial conditions" (vq^Vzo) at 
large z and integrating backwards. It is seen that the 
dependence on the particular choice of downflow condi- 
tions is very week since any phase trajectory quickly con- 
verges to the well defined stable manifold. Thus, even for 
a thread of moderate length the shape is uniquely deter- 
mined irrespective of the precise downstream conditions, 
just as we would expect. 

The asymptotic behaviour of the solution as z ^ cxd is 
easily seen to be controlled by only the two first terms in 
(EH), i.e. 



givmg 



i2z 



1 



h = .1- 

V 



■,1/2 



(2z) 



-1/4 



(24) 

(25) 
(26) 



This asymptotic solution is shown by the dot-dashed 
curve in Fig. ^ (marked " inertial" ) . 

The behaviour of the unstable manifold near the fixed 
point {v,Vz) = (0,0) can be found by expanding in z. 



Clearly v = Cz^ + 0{z-^) for the RHS of ^ to remain 
finite as z ^ 0. Inserting this expression into (|22|l . we 
see that the inertial term vvz can be neglected, since it 
contributes only as z^, whereas all other terms contribute 
with z°-terms, and we find 



with the (positive) solution 







C = 



1 + - yTT 



With this choice of C the solution 



(27) 
(28) 

(29) 
(30) 



is in fact an exact solution to l|22|l , when the inertial term 
vVz is neglected. This v{z) is shown by the dotted curve 
in Fig. n (marked "viscid"). 
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FIG. 1: The phase plane for the equation 1221 . It is seen that, 
upon backward integration, trajectories quickly converge to a 
well-defined "unstable manifold" (full line) for the fixed point 
(Wj^z) ~ (0,0). The asymptotic solution for large z, u ~ y^z, 
is shown dot-dashed and is governed by inertia and gravity. 
The asymptotic solution for small z, n ~ z^, is shown dotted 
and is obtained by neglecting inertia. 
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FIG. 2: Plot of the numerical solution of 1221 for different 
values of 7 



The crossover between the viscid and inertial solutions is 
roughly given by the value z* where they become equal, 
i.e. 



z* « const(^)2/3 



v2/3 



(32) 



On Fig. 121 we plot v{z) for various values of 7. 



IV. STABILITY PROPERTIES. 



For a very viscous fluid, where 75 ^ 1, the coefficient 



IS 



C 



1 

2^ 



(31) 



A. Stability of a fluid cylinder in the lubrication 
approximation 

Let us quickly go through the stability in this case 
using the dimensional equations H15|) - (|16|l instead of the 
dimensionless (|20ll - (|20ll since we want to take the limit 



4 



g = 0. In the absence of gravity term (the constant 1 on 
the RHS) where the stationary state is a cyhnder moving 
with constant velocity. This is the classical Rayleigh- 
Plateau instability in the long wave length approximation 
0- We thus assume 

V ^ vq + v{z,t), 
h — hQ + h{z,t). 

and obtain the linearized system 



ftp 



"•0 
1. 



(33) 



:hoVz 



It is convenient to go to the comoving frame 
y ^ z- vot, T = t. 

where 



hr 



1, 



iQVy. 



Transforming as usual to Fourier modes as 
{v, h) = (Ci, C2)exp{ixy + st), 
leads to the dispersion relation: 



1 



9 / 2ax'^ ^ , 
-vx^ ± \ — h v^x^ 



(34) 

(35) 
(36) 

(37) 
(38) 



which, within the long- wave region x <C /lo coincides with 
the well-known results for the classical Rayleigh-Plateau 
instability 1]. 

In the inviscid case ;/ = 0, we get (for the unstable 
mode with positive s): 



(39) 



B. Stability of the inviscid thread solution 

We now study directly the stability of the stationary 
states of H20|l - (|21|) in the limit of vanishing viscosity, i. e. 



Vt + VVz 



(40) 
(41) 



and we linearize around the stationary solution 
{vo{z),hoiz)) 1123-1211) as 



V = ■i;o(l + a) 
h = ho{l + b) 



(42) 



to obtain the linear system: 

at + voUz + 2voza = q^^^'^v^^^'^hz 



vobz 



vq 



+ \q-"W^^ozb (43) 

= (44) 



To get rid of the advection term, we introduce the 
stretched spatial variable y as 



y 



dz 



(45) 



so that VQ{z)dz — dy. We also define the function W{y) 
as 

W{y) = vo{z{y)) (46) 
and these definitions transform 1431) - (|44|l into 



at + ay + 2W ^Wya 



bt+by + ^ 



^q-^^^W-'/^Wyb (47) 
(48) 



For the inertial stationary solution vo{z) ~ \f2z we have 
explicitly 



^ = Y, my) = 2/ 



(49) 



and we finally transform ()47|l - (|48|l into the comoving 
frame of reference by 



to obtain 



y = X + t, t = t 



bt + ^ = 



(50) 
(51) 



at + 2a{x + t)-^ = q-^/^{x + t)-^^^b^ 



We now Fourier-transfom in x, assuming that the asymp- 
totic behaviour will not be influenced by the slow alge- 
braic variation with x in the denominators as long as 
i 3> a;, an assumption which will be verified in the Ap- 
pendix. Thus we find, in terms of the Fourier transforms 
d{k, t) and b{k, t), 

&t « -2at-i+ifcg-i/2i-3/2^+lq-i/2t-5/2^ (53) 



Making finally the substitution 

b^t-^B 



(54) 



5 



we find, retaining only the dominant term as t — > oo 
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The WKB ansatz 



gives 



B{k, t) = Bo exp( J $(fc, t')dt' 



B+{t) = eM2V2kq-^/H^/^) 



(55) 



(56) 



(57) 



Thus the typical instability time is tc q and using H23|) 
and I19() allows us to estimate the typical dimensional 
instability length - in time and in space, respectively - 



Zc 



gtl ^ ( 



-,2 -2 



Note that if h{z) 



and thus v{z) 



,2b 



(58) 
(59) 

the thread 



would be stable if2/7<6<l/2. This could be realized 
if the gravitational field increased as g{z) ~ z°' with a = 
46—1. The case studied above, h — 1/4, is slightly below 
the lower limit of stability. 



C. The effect of viscosity: spatially growing modes. 

With finite viscosity, linearization of (|2()(l and (|21(l in 
the same way as in the previous section equations leads 
to 

at + ay + 2W-^Wya ^ q-^'^W-^'^by 

+ W-^/'^Wyh 

2 

r-2 



+ {ayy + 2W-^Wyhy 

+ {Wy 0W-2W2 



2W-^W'^)a) (60) 



ht + by 







Perturbing again around the asymptotic state W 
is determined by H49|l . Again, we transform to the co- 
moving frame 150(1 . and in the regime t ^ x viscosity 
drops out - the viscid corrections are subdominant. But 
if we instead assume that x ^ t the situation is differ- 
ent. Thus we approximate by neglecting all explicit 
time-dependence (coming from W). It is thus natural 
to assume the following behaviour of the amplitudes of 
perturbation f[loj.pT|): 

(a(x, t), b{x, t)) = (A(x, s),B{x, s)) exp(st) (61) 

where s is a real number. This leads to the following 
system of equations: 



sB 



Ax 

~2 



(62) 



sA -t- 2W-^WxA 



5-1/2^-3/2^^ 



-f ^W~^ [Axx + 2W-^WA 



(63) 



After some manipulations we end up with the single sec- 
ond order equation: 



q 



2 2o 

- + -4 

X x^ 

-1/2^-1 



-1/2. 



(64) 



-3/2 



- JX 



"fX ^ 



Ax 



Let us now assume that we are looking for the solutions 
of (|64|l which grow in the positive direction of x (see 
[12]). From the physical point of view this means, that 
perturbations should remain finite in the area near the 
outlet. If the spatial coordinate satisfies the following 
condition: 



max(2s"\ 2i/S^/^s-^/4) < x < 47^52^ (65) 
the dominant terms in H64I) are give the simpler equation 



Ax 



g-i/27-is-i 



-x-'A = 



with the WKB-type solution 
oc exp 

which is valid for: 



max(7 ' s 



l/4e-l/4 ^-l/3«-l/3„-l/3^ 



-7 



< X 



(66) 



(67) 



(68) 



We conclude that viscosity gives rise to a superexponen- 
tial growth along the spatial variable x with a charac- 
teristic length li, ~ 7^^^^, which in dimensional variables 
becomes - Icj^^^ = a^l^'^ g-'^l^^v^l^ . 

For the large values of x the viscous effects drop out 
and leads to 



Axx + ^Ax + 2sq^^''x^l^A = (69) 

Ax 



which produces slowly decaying WKB-solutions: 

1/4-7/4^ (70) 



A^ cx x-i/4exp(±i\/2|s|gi/V/^) 



We can conclude, that in the first time instant t <C x the 
dynamics of the system is governed by spatially growing 
modes. But the region of validity for the "spatial asymp- 
totics" shrinks with time as 2: cx and the "temporal 
asymptotics" (cx exp(ti/^)) takes over. For large enough 
7, this regime will finally define the break-up of a fiow 
unless it have been destroyed already by spatially grow- 
ing modes (the latter seems to happen in the limit of 
small 7, e.g. for water). 
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V. DISCUSSION. 

The stationary flow of a long falling viscous thread, 
with the asymptotic shape h{z) ^ z'~^l'^ is unstable as ex- 
pected from the classical results for the stability of a fluid 
cylinder. The perturbations grow, however, very slowly, 
increasing asymptotically only as exp(const t"'^/^), where 
the constant is independant of viscosity. What is the 
role of viscosity then? First, with respect to stationary 
solution, viscosity defines the structure of the flow near 
the outlet. The crossover between viscous and inertial 
solutions is found to scale like oc ^/^Z"^. The less viscous 
threads are vulnerable due to spatial instability, since the 
perturbations grow like exp(j/~^/^a;^) along the thread. 
For large viscosities, the effects of spatial instability are 
weak, and the inviscid asymptotics will dominate the de- 
velopment of the break-up. It is interesting to note, that 
the final instability is so weak, that if the gravitational 
field was growing slowly i.e. g{z) ~ with 1/7 < a < 1, 
the thread would become asymptotically stable. 
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APPENDIX 

In this appendix we show that the explicit dependence 
on X in (|51|l can be neglected when i ^ Neglecting 
this variation led to the "local" solution (|57|l : 



ao = — oot 



(71) 
(72) 



where e — VSg and where we define the direct and 
inverse Fourier transforms as 
1 f 

a{k,t) — -z- I a{x,t) exp{—ikx)dx, (73) 



2tt 



a(x, t) 



a{k, t) exp{ikx)dk 



Let us now include effetcs of non-locality by Fourier 
transformation of (ISTl) 



1 f 2a exp(— ifca;) 
a* + — / : dx 



27r 



X + t 



_ _i/2j_ ( bx exp(-^fca:) 



1/2 1 f beKp(—ikx) 



dx 



(74) 



2Tr J 2(a;-ht)5/2 
Now we assume that |a;| t and expand integral ker- 
nels in (|7^ in power series of \x\/t. First we consider 
LHS of IZH): 



1 f 2a exp(—ikx) 
271 J x + t 



dx = 2t ^d 



2^i-"-iC!!i— / x"aexp(-ifca;)dx (75) 



We substitute: 



— / x"aexp{—ikx)dx = (i)" ( d 
2tt J \dk J 



(76) 



and get: 



1 f 2a exp(— ifcx) 



27r 



X + t 



dx — 2t ^d 



oo / J \ ^ 

2^t— c!^,(.)"y d 

n=l ^ ' 



(77) 



where C," is a binomial coefficient. Let us estimate cor- 
rection term in the RHS of H77|l for the local solution 
(|71(l - (|72|l . After some algebra we obtain: 



1 

2^ 



2ao(t) exp(ifca;) 



dx = 2r^do{t) 



x + t 

OO n 
n—1 m— 

X (-l)"-ifc-i-"+"e"t'^ (fee + (m - A)t-^''^) 



(78) 



Thus the highest power of t under the inner summa- 
tion is oc Taking into account general multi- 
plier oc t^"^^, in the leading order we have the correc- 
tion: oc t-(3n-i-7)/4^-igxp(j,^^i/4)^ This should be com- 
pared with the main contribution on the rhs. I|74fl . i.e. 
t-3/2f-iexp(efcii/4). Obviously: 



^-(3n+7)/4 ^ ^-3/2^ i OO, n=l,2. 



(79) 



Now we continue with the rhs. of 1)74(1 . Applying the 
same method, we get: 



1 

2^ 



bx exp(— i/cx) 
(x + ty/^ 



dx 



= ^ t-nC-,^,{^r+' ( ^ ) (kb) (80) 

For the correction term we get: 



oc t 



-3/2 



J2t-'^C-y,i^) 



n+l^n— 1 



ke + nt 



-1/4 



bo (81) 



Now the leading order term is oc i 3n/4^ 3/2^^j^^-j^ which 
is dominated by the main term oc t^^^'^bo{t) when t —> oo. 
As for the last term in the RHS of H74() . it only produces 
minor corrections to the main solution even in the leading 
order of magnitude. Thus we see that for \x\/t ^ 1 we 
can completely neglect the effects of non- locality in (|74|) . 
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